clc
clear
currentDir = pwd; % 获取当前目录路径
parentDir = fileparts(currentDir); % 获取当前目录的上一级目录路径
% grandParentDir = fileparts(parentDir); % 获取当前目录的上两级目录路径
% greatGrandParentDir = fileparts(grandParentDir); % 获取当前目录的上三级目录路径

%指定需要查找的文件夹
subdir = 'DATA';
%指定需要查找的文件名
filename = '413168000(NoRepeated).mat';%413168000.mat'; %matlab中包含的是剔除
%构造需要查找的文件的完整路径
full_path = fullfile(parentDir, subdir, filename);

load( full_path)  %此处使用的数据是原始org数据经过后面的[]处理后得到的数据

%% 数据预处理
% AIS(AIS.Speed<1,:)=[];%锚泊状态不处理
AISData=data;

% 设置参数
MinTime = 30; % 时间停留门限，单位为秒
MaxSpd = 0.5; % 速度停留门限，单位为km/h
MinPts = 2; % 邻域内的最小轨迹点数

% 数据预处理 - 剔除停留点
processed_data = AISData;
stay_points = []; % 初始化停留点数组
for i = 1:height(AISData)
    if i == 1 || (numel(stay_points) == 0 || stay_points(end) ~= i-1 && abs(AISData.Speed(i)) < MaxSpd && abs(AISData.Speed(i-1)) < MaxSpd)
        % 检查是否为停留点
        stay_time = 0;
        while i < height(AISData) && abs(AISData.Speed(i)) < MaxSpd && abs(AISData.Speed(i-1)) < MaxSpd
            stay_time = stay_time + seconds(AISData.Time(i) - AISData.Time(i-1));
            i = i + 1;
        end
        if stay_time > MinTime
            stay_points = [stay_points, i-1]; % 标记停留点
        end
    end
end

% 剔除停留点
processed_data(stay_points, :) = [];

% 异常点识别 - 网格划分
% 定义网格大小
grid_size = 0.01; % 网格大小，单位为度

% 将经纬度转换为网格索引
lon_idx = floor(processed_data.Lon / grid_size);
lat_idx = floor(processed_data.Lat / grid_size);

% 计算每个网格的点数
[unique_lon, ~, lon_idx] = unique(lon_idx);
[unique_lat, ~, lat_idx] = unique(lat_idx);
[grid_coords, ~, grid_idx] = unique([lon_idx lat_idx], 'rows');
grid_count = accumarray(grid_idx, 1);

% 定义异常点阈值
outlier_threshold = 8; % 可以根据实际情况调整

% 识别异常点
outlier_grids = find(grid_count < outlier_threshold);
outliers_idx = false(height(processed_data), 1);
for k = 1:length(outlier_grids)
    grid_idx = outlier_grids(k);
    lon_grid = grid_coords(grid_idx, 1);
    lat_grid = grid_coords(grid_idx, 2);
    points_in_grid = (lon_idx == lon_grid) & (lat_idx == lat_grid);
    outliers_idx(points_in_grid) = true;
end


% 剔除异常点
cleanedAISData = processed_data(~outliers_idx, :);

% 显示结果
disp('剔除异常点后的轨迹数据：');
disp(cleanedAISData);

% 可视化结果
figure;
plot(processed_data.Lon, processed_data.Lat, 'b.'); % 原始轨迹
hold on;
plot(cleanedAISData.Lon, cleanedAISData.Lat, 'r.'); % 剔除异常点后的轨迹
legend('原始轨迹', '剔除异常点后的轨迹');
xlabel('Longitude');
ylabel('Latitude');
title('AIS Trajectory Anomaly Detection');
grid on;


%计算偏转角标准差,
lhm.degrees_org = calculate_track_degree(table2array(AISData(:,{'Lon','Lat'})));
lhm.degrees_deleted = calculate_track_degree(table2array(cleanedAISData(:,{'Lon','Lat'})));

lhm.mean_org=mean(lhm.degrees_org);
lhm.mean_deleted=mean(lhm.degrees_deleted);
lhm.std_org=std(lhm.degrees_org);
lhm.std_deleted = std(lhm.degrees_deleted);